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STABILITY  MEASURES  FOR  SPECTRAL  ANALYSIS  USING 
DISCRETE  SAMPLING  WITH  THE  KAISER-BESSEL 
OR  DOLPH-CHEBYSHEV  WINDOW 


1.  INTRODUCTION 

Power  spectral  estimation  from  a  sample  of  a  stationary  random  process 
has  important  detection  and  estimation  applications  in  sonar.  Consequently, 
we  would  like  to  obtain  as  accurate  an  estimate  as  possible  (accurate  meaning 
that  the  estimator  has  a  low  standard  deviation  relative  to  its  mean) .  This 
measure  of  estimator  reliability  is  known  as  stability,  ideally  then,  we  would 
like  an  estimate  that  is  as  stable  as  possible,  given  a  finite  sample  of  data. 

The  power  spectrum  estimator  examined  in  this  study  is  the  standard 
overlapped  fast  Fourier  transform  (FFT)  processor.  With  this  method,  ashd- 
ing  weight  function  is  applied  multipUcatively  to  the  data.  The  magnitude- 
squared  FFTs  of  each  weighted  data  segment  axe  then  summed  to  form  the 
overall  estimate. 


The  Kaiser-Bessel  and  Dolph-Chebyshev  functions  axe  used  for  the  shd- 
ing  weight.  The  Kaiser-Bessel  window  is  characterized  by  a  free  parame¬ 
ter  13  that  controls  the  peak  sidelobe  level  and  mainlobe  width;  the  Dolph- 
Chebyshev  window  has  the  narrowest  main  lobe  for  a  given  peak  sidelobe 
level.  Both  windows  are  thoroughly  discussed  in  the  classic  Harris  paper.  ^ 

Our  primary  objective  is  to  determine  the  most  stable  estimate  that  can 
be  obtained  with  the  two  different  weightings;  that  is,  what  values  of  overlap 
(85,  90,  98  percent,  etc.)  will  yield  the  highest  stability  measure?  Rules  of 
thxnnb  are  developed  for  establishing  the  optimal  amount  of  overlap,  based 
upon  the  two  windows  and  their  peak  sidelobe  levels.  Moreover,  since  a  50- 
percent  overlap  is  often  used  in  practice,  we  derive  the  exact  performance 
of  this  particular  choice  of  processing  (in  terms  of  stabihty)  relative  to  the 
optimum  percentage  overlap. 

At  the  conclusion  of  this  investigation,  we  find  that  although  a  50- 
percent  overlap  is  usually  suboptimal,  it  often  yields  near-optimal  stability. 
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2.  DEFINITION  AND  INTERPRETATION  OF 
EQUIVALENT  DEGREES  OF  FREEDOM 


Suppose  that  we  have  a  set  of  J  statistically  independent,  zero-mean  Gaus¬ 
sian  random  variables  of  common  unknown  power  which  we  wish  to  esti¬ 
mate.  The  best  that  we  can  do  in  this  case  is  to  compute  the  power  estimate 


(1) 


This  quantity  <Tg  is  an  unbiased  estimate  of  since  its  mean,  Av{(Tg},  is  obviously 
<7^,  the  desired  quantity.  At  the  same  time,  since  variances  of  independent  random 
variables  add,  the  variance  of  the  power  estimate  is 

Var{<T2}  =  jVar{p|}.  (2) 

And,  since  the  fourth  moment  of  Gaussian  random  variable  is  the  quantity 
in  equation  (2)  is  simply 

Var{<Tgn  =  ^  ,  (3) 

which  decays  with  J,  the  number  of  independent  samples  available.  Naturally,  we 
want  J  to  be  as  large  as  possible,  in  order  that  the  variance  of  the  estimate  <Tg  be 
at  a  minimum. 


Now  consider  the  quality  measure  of  estimate  <Tg,  defined  as  the  equivalent 
degrees  of  freedom,^  namely. 


EDF  = 


2Av2{<t2} 

Var{o-2} 


r 

2(7V  J  ’ 


(4) 


where  we  have  used  equation  (3).  This  ratio  gives  exactly  the  number  of  statistically 
independent  samples  available  in  the  original  data  set  {gf^}.  Therefore,  it  is  natural 
and  useful  to  consider  extending  the  EDF  definition  in  equation  (4)  to  an  arbitrary 
random  variable  v  (formed  perhaps  as  a  sum  of  partially  dependent  random  vari¬ 
ables)  as  a  measure  of  the  effective  number  of  independent  contributors  to  v.  That 
is,  we  define,  in  general, 


EDF  = 


2Av"{t;} 

Var{v} 


(5) 
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for  random  variable  v.  This  statistic  is  attractive  in  that  it  only  requires  the  cal¬ 
culation  of  the  two  lowest  moments  of  v.  (If  the  mean  of  v  is  zero,  the  definition 
in  equation  (5)  will  have  to  be  discarded;  it  is  most  meaningful  when  applied  to 
estimates  of  power  quantities,  which  are  inherently  positive.) 

If  random  variable  v  depends  on  some  parameter  X,  a  search  can  be  con¬ 
ducted  on  that  parameter  in  an  effort  to  maximize  EDF,  that  is,  minimize  the 
standard  deviation  of  v  relative  to  its  mean.  In  this  way,  the  relative  stability  of 
estimate  v  can  be  optimized,  subject,  of  course,  to  whatever  constraints  might  have 
to  be  imposed  on  parameter  X  due  to  practical  limitations. 
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3.  STABILITY  OF  OVERLAP  SPECTRAL 
ANALYSIS  FOR  DISCRETE  SAMPLING 


Stationary  random  process  x{t)  has  correlation  function  Rx{r)  =  x{t)x{t  —  r) 
and  (double-sided)  power  density  spectrum  Gx(f).  This  process  x{t)  is  sampled  at 
increments  of  A  seconds,  where  A  should  be  chosen  small  enough  to  keep  spectral 
abasing  under  control;  if  not,  we  shall  be  estimating  an  excessively  aliased  spectrum 
instead  of  Gx{f)*  The  data  samples  available  are  x(0),  a;(A), ...  ,x((T  —  1)A),  for 
a  total  of  T  samples. 

The  total  data  record  of  T  samples  is  broken  up  into  segments  of  size  N, 
which  need  not  be  a  power  of  2.  Furthermore,  these  segments  may  overlap  in  time, 
thereby  resulting  in  some  data  samples  being  processed  in  multiple  segments.  In 
addition,  each  data  segment  is  weighted  multiplicatively  with  a  set  of  (delayed) 
coefficients,  w(0),tt;(l), . . .  ,w{N  —  1),  prior  to  being  Fourier  transformed,  in  order 
to  control  the  sidelobe  behavior  of  the  windowed  Fourier  transform  W  (/)  associated 
with  the  set  of  weights.  That  is,  window 


N-l 

W (/)  =  ^  w{n)  exp(— z27r/An)  for  all  /  ,  (6) 

n=0 

which  has  period  1  /A  in  frequency  /. 

The  segmenting  and  weighting  procedure  applied  to  the  available  data  is 
illustrated  in  figures  1  and  2.  The  total  number  of  pieces  averaged  in  obtaining  the 
power  density  spectrum  estimate  is  P,  which  can  vary  over  a  wide  range  (under 
our  control).  The  pth  piece  of  weighted  data  is  subjected  to  the  (integer)  shift  Sp, 
where  1  <  p  <  P.  In  particular,  we  will  choose  =  0,  in  order  to  cover  the  first 
data  sample  a;(0).  Also,  we  will  choose  Sp  =  T  -  N,  so  that  the  last  data  segment 
covers  the  last  available  data  sample  x{{T  -  1)A).  Figure  1  shows  the  case  where 
the  data  segments  do  not  overlap;  figure  2  shows  the  case  where  the  data  segments 
do  overlap  because  a  large  number  of  segments  (large  P)  have  been  used.  Note 
also  from  figure  2  that  all  the  data  points  may  not  be  used,  depending  upon  the 
relationship  between  N,  T,  and  P. 

‘This  requirement  for  the  sampling  rate  A  is  equivalent  to  the  Nyquist  sampling  requirement 
in  the  frequency  domain. 
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Figure  1.  Weighted  and  Overlapped  Data  Segments  (Nonoverlapping  Case) 


large  P 


12  3  P 


Figure  2.  Weighted  and  Overlapped  Data  Segments  (Overlapping  Case) 


In  order  to  maximize  the  stability  (minimize  the  variance)  of  the  spectral 
estimate,  we  want  to  separate  the  P  individual  data  segments  as  much  as  possible; 
that  is,  we  want  to  minimize  the  amount  of  overlap  of  the  individual  segments.  We 
accomplish  this  by  first  defining  the  (noninteger)  quantity 

p  - 1  ’ 

multiples  of  this  shift  would  result  in  aU  the  P  segments  being  equally  spaced. 
However,  since  Sc  and  its  multiples  are  not  integers,  we  will  adopt  the  closest  integer 
spacing  for  the  shifts  {Sp}  according  to  the  rule 

Sp=  \{p-  l)^^  for  1  <  p  <  P ,  (8) 

where  [t?]  denotes  the  closest  rounded  integer  to  v.  As  special  cases,  we  have  Pi  =  0 
andSp  =  T-N. 

For  arbitrary  analysis  frequency  /fc,  the  pth  voltage  density  estimate  (using 
weighted  sampled  overlapped  data)  is  defined  as  the  complex  quantity 

Vpifk)  =  ^  x{nA)  w{n-  Sp)  exp  (-i27r/fcnA)  for  1  <  p  <  P ;  (9) 
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the  limited  nonzero  range  of  the  weights,  w(0)  to  w{N  -  1),  will  automatically 
terminate  the  infinite  summation  in  equation  (9)  at  limits  Sp  to  Sp  +  N  —  1,  as 
shown  in  figure  1.  (If  we  choose  to  take  fk  =  k/{KA),  where  K  is  a  power  of  2, 
then  equation  (9)  can  be  accompHshed  as  a  iiT-point  FFT.  Also,  FFT  size  K  can  be 
chosen  larger  than  segment  length  to  realize  interpolated  values  of  the  spectral 
estimate  at  frequency  increment  {KA)~^  of  specified  fineness.) 

The  power  density  spectrum  estimate  at  frequency  fk  is  now  derived  from 
equation  (9)  as  the  finite  average 

=  (10) 

p=i 


We  are  interested  in  the  mean  and  standard  deviation  of  random  variable  Y{fk) 
and,  in  particular,  in  the  stability  measure  known  as  the  EDF,  namely. 


EDF  = 


2Av^{y(/fe)} 

Var{r(A)}  • 


(11) 


For  given  values  of  T  and  N,  the  value  of  EDF  will  depend  on  P,  the  number 
of  pieces  in  figure  1  and  in  the  estimate  of  equation  (10).  We  wish  to  make  EDF 
as  large  as  possible  to  minimize  the  variance  and  yet  maintain  P  within  reasonable 
bounds.  Too  small  a  value  of  P  results  in  incomplete  use  of  the  available  data;  see 
figure  1  for  an  illustration  of  the  nonoverlapping  case.  On  the  other  hand,  too  large 
a  value  of  P  results  in  excessive  overlap  and  an  attendant  amount  of  unnecessary 
calculations  with  little  additional  reduction  in  variance.  The  average  fractional 
overlap  (FO)  of  the  segments  is  available  from  equation  (7)  as 

FO  =  max(:^^^,o)=max(l-|,o).  (12) 

For  example,  assume  that  T  =  1024  total  data  points  and  that  the  data 
segments  consist  of  =  128  data  points.  If  P  =  8  segments  were  desired,  then 
Sc  =  (1024  —  128)/(8  —  1)  =  896/7  =  128  from  equation  (7)  and  FO  =  0  from 
equation  (12).  This  result  clearly  shows  that  P  =  8  segments  of  128  points  each 

^  The  additional  K—N  data  points  all  have  zero  value;  this  procedure  is  known  as  zero  padding. 
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can  be  obtained  from  1024  total  points  without  any  overlap.  If  P  is  now  changed 
from  8  to  16,  then  Sc  becomes  896/15  59.67,  and  FO  becomes  0.53,  since  with 

16  pieces  of  128  points  each,  there  must  be  overlap  if  there  are  only  1024  total  data 
points  available.* 


The  EDF  is  evaluated  in  appendix  A,  with  the  result  that 

2P2 


EDF  = 


E 

p,q=l 


(13) 


where  (pwi’iTT-),  the  autocorrelation  of  the  weight  sequence,  is 
(py;{m)  =  w{n)  w{n  —  m)  for  all  m . 


(14) 


This  sequence,  {^^(m)},  is  zero  for  |m|  >  N,  since  {w{n)}  is  of  length  N.  Shifts 
{S'p}  are  available  from  equations  (7)  and  (8)  for  each  choice  of  P. 


The  detailed  spectrum  Gx{f)  of  random  process  x{t)  does  not  appear  in  the 
end  result  for  the  EDF  in  equation  (13),  nor  does  the  particular  analysis  frequency 
fk-  The  reason  for  this  behavior  is  our  presumption  that  the  magnitude-squared 
window  |W(/)|^  has  a  mainlobe  width  which  is  narrower  than  the  detail  in  the 
spectrum  Gxif)  in  the  neighborhood  of  fk  and  that  the  sidelobes  of  the  window 
W{f)  in  equation  (6)  are  rather  small. ^  Both  these  properties  can  be  achieved  by 
the  selection  of  a  large  enough  segment  duration,  NA,  and  by  the  use  of  windows, 
W (/),  with  good  sidelobes,  such  as  can  be  achieved  from  the  Kaiser-Bessel  or  Dolph- 
Chebyshev  weightings  {ty(n)}. 


We  can  now  concentrate  on  maximizing  the  quality  ratio  EDF  in  equa¬ 
tion  (13)  by  the  choice  of  the  number  of  pieces  P,  knowing  that  this  will  minimize 
the  standard  deviation  of  the  spectral  estimate  Y{fk)  relative  to  its  mean  value. 

must  be  noted  that  the  amount  of  FO  available  in  actual  sonar  systems  is  often  restricted 
to  values  such  as  0.5, 0.75,  and  0.875.  Thus,  an  FO  of  0.53,  as  discussed  in  the  last  example,  may  not  be 
possible  in  practice. 

^This  is  a  reasonable  assumption  for  practical  FFT  spectral  estimation  processes. 
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4.  RULES  OF  THUMB 


The  following  tables  show  the  optimum  overlap  for  maximizing  stability  mea¬ 
sure  EDF  and  indicate  what  fraction  of  this  maximal  amount  is  achieved  for  50- 
percent  overlap.  From  these  tables,  it  can  be  seen  that  the  following  rules  of  thumb 
apply. 


First,  tables  1  and  2  show  in  the  fifth  column  that  the  optimum  overlap 
required  to  achieve  maximum  EDF  is  directly  affected  by  the  desired  sidelobe  level, 
data  length  T,  and  segment  size  N.  As  the  sidelobe  level  decreases,  the  amount  of 
overlap  required  for  maximum  EDF  increases.  These  two  tables  also  indicate  (see 
the  eighth  column)  that  for  both  the  Dolph-Chebyshev  and  Kaiser-Bessel  windows, 
50-percent  overlap  provides  near  optimum  EDF,  especially  for  shallower  sidelobe 
levels  (selectable  by  the  parameter  /?  for  Kaiser-Bessel). 

Figures  3  through  6  also  illustrate  a  rule  of  thumb,  which  relates  to  how 
the  EDF  changes  as  a  function  of  overlap  for  the  Kaiser-Bessel  weightings  and  four 
different  peak  sidelobe  levels.  For  example,  figure  4  shows  EDF  as  a  function  of  the 
number  of  pieces  (sections)  for  a  Kaiser-Bessel  window  designed  for  a  -40  dB  peak 
sidelobe  level  with  T  =  1024  and  N  =  128.  For  this  example,  the  maximum  EDF 
value  of  29.42  is  attained  at  FO  =  0.632  (20  sections).  The  EDF  value  at  50-percent 
overlap  is  27.82,  resulting  in  a  loss  of  only  5.4  percent  relative  to  the  maximum  EDF 
value  obtained  at  the  optimal  overlap.  It  is  obvious  that  once  the  knee  of  the  curve 
has  been  attained,  further  increases  in  the  fractional  overlap  have  little  effect  on 
improving  EDF  (i.e.,  obtaining  better  spectral  estimates). 

Another  general  rule  can  be  summarized  as  follows:  for  a  given  N  and  T, 
EDF  increases  with  deeper  sidelobe  levels.  Thus,  we  expect  a  larger  value  of  EDF 
for  a  peak  sidelobe  level  of  -60  dB  than  for  a  peak  sidelobe  level  of  -30  dB.  The 
vahdity  of  this  last  statement  can  be  verified  by  examination  of  tables  1  and  2. 

It  is  helpful  to  have  a  simple  rule  of  thumb  to  obtain  the  maximum  value 
of  EDF  that  can  be  realized  for  the  two  windows  considered  here.  Accordingly,  we 
have  extracted  the  following  approximation  from  the  examples  in  tables  1  and  2. 
It  is  based  upon  reference  3  and  is  a  function  of  the  ratio  T/N,  namely,  the  ratio 
of  the  total  number  of  data  points  to  the  length  of  the  segments  employed.  It  was 
obtained  by  fitting  the  empirical  examples  in  tables  1  and  2  by  a  straight  line  in  the 
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Table  1.  Optimum  Overlap,  Maximum  EDF,  and  EDF  Attainable  at 
50-Percent  Overlap  for  Kaiser-Bessel  Weighting 


SLL 

(dB) 

T 

N 

T/N 

Optimum 

Overlap 

Max 

EDF 

Max  EDF 
T/N 

Fraction  at 
50%  Overlap 

-30 

512 

128 

4 

0.57 

12.2 

3.05 

0.99 

8 

0.59 

mJmM 

3.26 

0.99 

wmim 

16 

0.61 

wy;i 

3.36 

0.98 

-30 

2048 

256 

8 

0.59 

26.0 

3.25 

0.99 

-30 

128 

32 

0.62 

109.2 

3.41 

0.98 

-30 

1024 

8 

0.59 

25.9 

0.99 

-40 

512 

128 

4 

0.63 

BU 

3.42 

0.96 

-40 

1024 

128 

8 

0.63 

0.95 

2048 

128 

16 

0.66 

61.0 

3.81 

0.94 

2048 

1^1 

8 

0.63 

29.3 

3.67 

0.95 

-40 

4096 

32 

0.67 

124.1 

3.88 

0.94 

-40 

8192 

1024 

8 

0.63 

29.3 

3.66 

0.95 

-50 

512 

128 

4 

0.67 

14.9 

3.73 

0.91 

-50 

1024 

128 

8 

0.68 

32.3 

4.04 

0.89 

2048 

16 

bq^hi 

BEISBi 

0.89 

2048 

IBS 

8 

bb^b 

-50 

4096 

128 

32 

0.71 

137.0 

4.28 

0.88 

jggl 

8192 

1024 

8 

32.1 

512 

128 

4 

IBBHI 

0.86 

-60 

1024 

128 

8 

0.71 

34.9 

4.37 

0.84 

-60 

2048 

128 

16 

0.73 

72.8 

4.55 

0.84 

-60 

2048 

256 

8 

0.71 

34.8 

4.35 

.  . 

0.85 

-60 

4096 

0.74 

IBS 

4.64 

0.83 

IfiTigl 

8 

0.71 

BU 

4.34 

0.85 
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Table  2.  Optimum  Overlap,  Maximum  EDF,  and  EDF  Attainable  at 
50-Percent  Overlap  for  Dolph-Chebyshev  Weighting 


SLL 

(dB) 

T 

N 

T/N 

Optimum 

Overlap 

Max 

EDF 

Max  EDF 
T/N 

Fraction  at 
50%  Overlap 

-30 

512 

128 

4 

0.57 

12.2 

3.1 

0.99 

-40 

1024 

128 

8 

0.59 

26.4 

3.3 

0.98 

-40 

■afeVJ 

1024 

8 

3.4 

0.96 

-40 

128 

16 

0.63 

3.4 

0.98 

-40 

16384 

1024 

■a 

0.84 

54.5 

3.4 

0.95 

-40 

4096 

128 

0.65 

110.6 

0.98 

-40 

32768 

32 

0.74 

114.0 

3.6 

0.95 

8 

0.63 

HQ 

3.6 

mm 

16 

0.66 

BitlKl 

3.8 

-50 

4096 

0.67 

122.7 

3.8 

0.94 

-50 

2048 

mm 

8 

0.63 

29.0 

0.95 

-50 

2048 

4 

0.63 

13.5 

3.4 

0.96 

-50 

mm 

8 

0.63 

29.0 

3.6 

0.95 

mm 

0.66 

60.3 

3.8 

0.95 

«!■ 

miozM 

16 

0.66 

60.0 

0.95 

32 

0.67 

122.7 

3.8 

0.94 

-50 

32768 

1024 

32 

0.68 

122.1 

3.8 

0.94 

128 

8 

0.91 

-60 

1024 

8 

31.4 

3.9 

0.91 

-60 

2048 

128 

16 

65.7 

4.1 

0.90 

-60 

16384 

1024 

16 

65.3 

4.1 

0.90 

-60 

4096 

128 

32 

0.70 

133.9 

4.2 

0.90 

-60 

32768 

mm 

0.70 

133.1 

4.2 

-70 

1024 

■wa 

8 

4.3 

-70 

■acfei 

■riwi 

8 

33.8 

4.2 

0.87 

mm 

16 

lIBBHi 

70.7 

4.4 

0.86 

-70 

16384 

1024 

16 

0.71 

70.3 

4.4 

0.86 

-70 

4096 

128 

32 

144.2 

4.5 

0.85 

-70 

32768 

mum 

■KUH 

143.4 

4.5 

0.86 
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Figure  4.  Kaiser-Bessel  Window  With  Peak  Sidelobe  Level  of  -fO  dB 
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EDF 


Percsentage  overlap 


Figure  5.  Kaiser-Bessel  Window  With  Peak  Sidelohe  Level  of  ~50  dB 


Figure  6.  Kaiser-Bessel  Window  With  Peak  Sidelohe  Level  of  -60  dB 


ratio  T/N.  We  find  that 

max  EDF  ^  ^ .  (15) 

where  77  and  6  depend  on  the  type  of  window  and  the  decibel  sidelobe  level  according 
to  tables  3  and  4.  This  rather  accurate  approximation  has  been  verified  by  examples. 


Table  3.  Values  ofr]  and  0  in  the  Approximation  of  Maximum  EDF  by 
(^  —  0)  77  for  the  Kaiser-Bessel  Weights 


Sidelobe 
Level  (dB) 

Values  of 

V 

Values  of 

e 

-40 

3.93 

0.54 

-50 

4.35 

0.59 

-60 

4.72 

0.63 

-70 

5.07 

0.67 

Table  4.  Values  of  77  and  9  in  the  Approximation  of  Maximum  EDF  by 
(^  —  ^)  77  for  the  Dolph-Chebyshev  Weights 


Sidelobe 
Level  (dB) 

Values  of 

V 

Values  of 

9 

-40 

3.51 

0.49 

-50 

3.89 

0.54 

-60 

4.25 

0.58 

-70 

4.58 

0.61 
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5.  SUMMARY 


In  this  report,  we  have  investigated  how  to  determine  the  optimal  amount 
of  data  overlap  to  produce  spectral  estimates  of  maximum  stability.  For  the  two 
types  of  spectral  windows  examined,  namely,  Kaiser-Bessel  and  Dolph-Chebyshev, 
the  optimal  amount  of  overlap  varied  from  59  to  93  percent  and  depended  greatly  on 
the  ratio  T/N,  the  ratio  of  the  total  number  of  data  points  available  to  the  number 
of  data  points  in  a  segment.  In  spite  of  this  dependency,  it  was  shown  through 
extensive  computation  that  a  50-percent  overlap  results  in  a  loss  of  stability  (as 
measured  by  the  EDF  value)  of  only  1  to  15  percent,  with  the  loss  frequently  less 
than  10  percent. 
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APPENDIX  A 

DERIVATION  OF  MEAN  AND  VARIANCE 

The  power  density  spectrum  estimate  is  given  by  equation  (10)  in  conjunc¬ 
tion  with  equation  (9).  For  later  use,  the  crosscorrelation  of  two  different  spectral 
estimates  at  general  analysis  frequency  fk  is  given  by  the  complex  ensemble  average 

Cpq  =  yp{fk)y*g{fk) 

=  ^  ^  x{nA)x*{mA)  w{n  —  Sp)  w*{m  —  Sg)  exp  l—i27r fiiA{n  —  m)]  , 

n  m 

(A-l) 

where  we  have  allowed  for  complex  processes  x{t)  and  complex  weights  {w(n)}. 
When  we  express  correlation  Rx{nA—mA)  in  equation  (A-l)  in  terms  of  its  spectrum 
Gx{f)  according  to  a  Fourier  transform,  the  average  becomes 

Cpg  =  df  Gx{f)  ^  exp  (i27r/An)  w{n  —  Sp)  exp  (— i27r/feAn) 

n 

X  ^^exp  (— z27r/Am)  w*{m  —  Sg)  exp  {i2-KfkAm) 

m 

=  jdf  G.{/)  \W(ft  -  ff  exp  [j2,r(S,  -  Sp)A(A  -  /)]  ,  (A-2) 

where  we  have  used  the  window  definition  from  equation  (6).  Note  that  the  only 
relevant  feature  of  the  window  is  its  magnitude-squared  value,  \W{f)f,  insofar  as 
the  crosscorrelation  Cpg  is  concerned. 

All  the  functions  in  equation  (A-2)  have  period  1/A  in  /,  except  for  spectrum 
Gx(,f),  which  decays  with  /  (and  leads  to  convergence  of  the  integral).  Therefore, 
we  can  express  crosscorrelation  Cpg  in  terms  of  the  (double-sided)  aliased  spectrum 

G,(/)  =  ^G,(/-^)  for  all/  (A-3) 

m 

(which  also  has  period  1/A  in  /)  according  to 

C„=  I  df  Gaif)  mh  -  /)p  exp  (i2,r(S,  -  S,)A(/s  -  /)]  ,  (A-4) 

1/A 


A-l 


where  we  can  now  integrate  over  any  period  of  length  1/A.  This  is  an  exact  result. 
Thus,  it  is  always  the  aliased  spectrum  Ga{f)  that  governs  the  statistics  of  the 
spectral  estimate  Y{fk)  in  equation  (10).  This  is  a  natural  consequence  of  process 
x{t)  being  sampled  at  time  increment  A,  and  is  not  due  to  the  particular  method 
of  spectral  analysis  employed.* 

If  the  original  stationary  random  process  a3(t)  is  prefiltered  by  an  analogue 
filter  prior  to  being  sampled  at  time  increment  A  seconds,  the  frequency  components 
in  Gx{f)  outside  the  range  —0.5/A  to  0.5/A  Hz  can  be  sufficiently  suppressed  so  that 
the  effects  of  aliasing  are  negligible.  Thus,  although  there  are  always  replications  of 
Gx{f)  in  Gaif)  at  spacing  1/A  Hz,  through  proper  prefiltering,  their  overlap  will 
be  insignificant  and  abasing  effects  can  be  ignored. 

A  special  case  of  equation  (A-4)  is  given  by  the  autocorrelation  value 
Cn.=  j  dfG.(f)  mh  -  f)f  ,  {A-5) 

1/A 

which  is  recognized  as  the  convolution  (at  analysis  frequency  /fc)  of  the  aliased 
spectrum  with  the  magnitude-squared  window,  the  integral  being  conducted  over 
any  period  of  length  1/A. 

At  this  point,  we  assume  that  the  magnitude-squared  window  \W{f)\^  is 
narrow  enough  so  that  only  the  value  of  the  aliased  spectrum  Ga{f)  at  analysis 
frequency  fk  matters  in  the  integral  of  equation  (A-4).  That  is,  the  detail  in  Ga{f) 
near  fk  is  broader  than  the  window  width.  Then,  we  obtain  the  good  approximation 

c„  =  G.iA)  j  df  \W{h  -  /)p  exp  [i27r(S,  -  S',) A(A  -  /)]  .  (A-6) 

1/A 

The  integral  in  this  equation  can  be  simplified  by  the  substitution  of  the  window 
definition  from  equation  (6)  and  the  interchange  of  the  integral  with  the  two  sum¬ 
mations,  yielding  the  compact  result 

Cp,  =  i  Ga{fk)  K  {S,  -  5p) ,  (A-7) 

^Sampling  x{t)  at  time  increment  A  causes  its  power  spectrum  to  be  periodically  reproduced 
with  a  period  of  1/A.  If  A  is  not  small  enough,  these  repetitions  will  overlap,  causing  aliasing. 


A-2 


where  is  the  autocorrelation  of  the  weight  sequence: 

w{n)  w* (n  —  m)  for  all  m .  (A-8) 

n 

In  particular,  equation  (A-7)  yields  Cpp  =  ^  Gaifk)  (f>w{0)  for  1  <  p  <  P. 

We  can  now  address  the  statistics  of  power  spectral  estimate  Y (fk)  in  equa¬ 
tion  (10).  By  use  of  equation  (A-1),  we  have  the  mean  value 

1  ^  1 

Av{y(A)}  =  p  =  -  G.{fk)M0)  ■  (a-9) 

P=1 

This  result  holds  regardless  of  the  statistics  of  process  x{t);  that  is,  x{t)  need  not  be 
Gaussian.  Since  the  absolute  scaling  is  arbitrary,  we  can  select  <^^,(0)  =  A  so  that 
the  mean  value  equals  the  desired  quantity,  namely,  spectral  value  Ga{fk)-  Then, 
estimate  Y{fk)  would  be  unbiased. 

To  evaluate  the  variance  of  Y (/t),  we  need  to  assume  that  the  random  vari¬ 
ables  {Vpifk)}  in  equation  (9)  are  zero-mean  Gaussian;  this  would  happen,  for 
example,  if  input  process  x{t)  were  zero-mean  Gaussian.  The  mean  square  value  of 
Y {fk)  is  then  given  by 

_  I  P  _ 

=  -p^Y^  ypylyqVq 

p,g=l 

1  P _ _ _ 

^  S  (ypy*p  y^y*<i + ypy<i  y^yq + ypy^  ypyq)  > 

p,q=l 

where  we  have  used  the  Gaussian  character  of  {y^f.  The  first  term  in  equation 
(A-10)  leads  to  the  square  of  the  mean  value  in  equation  (A-9).  The  second  term 
in  equation  (A-10)  is  approximately  zero,  as  may  be  confirmed  by  conducting  an 
analysis  similar  to  that  of  equation  (A-2)  for  instead  of  a  magnitude-squared 
window,  we  obtain  two  displaced  windows,  W{fk  —  f)W{fk  +  /),  which  are  widely 
spaced  in  frequency  (except  when  fk  is  equal  to  a  multiple  of  the  Nyquist  frequency 
^),  thereby  resulting  in  a  small  output.  The  third  term  in  equation  (A-10)  therefore 
leads  to  the  variance  in  the  form 

Vs^{Y{M]  =  \c„\^  =  (S,  -  ,  (A-11) 

p,?=l  p,g=i 


A-3 


where  we  have  used  equation  (A-7). 


The  EDF,  defined  in  equation  (11),  is  now  available  from  equations  (A-9) 
and  (A-11)  as 


p 

E 


2P^ 

(t>wiSq  -  Sp)  ^ 

0W)(O) 


(A-12) 


where  the  exact  values  of  the  aliased  spectrum  Gaifk),  the  particular  analysis  fre¬ 
quency  fk,  and  the  samphng  increment  A  have  canceled  out.  Thus,  the  relative 
stability  of  the  power  density  spectrum  estimate  Y{fk)  (of  the  aliased  spectrum) 
depends  only  on  the  shape  (not  the  level)  of  the  weight  sequence  {w{n)}  and  integer 
parameter  values  N,  T,  and  P.  The  P  integer  shifts  {Sp}  are  given  by  equations  (7) 
and  (8),  and  depend  on  the  particular  number  of  pieces,  P,  being  considered;  see 
figures  1  and  2  for  the  method  of  data  processing  used  here. 
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APPENDIX  B 

SOME  BANDWIDTH  MEASURES  FOR  DISCRETE  WEIGHTING 


Various  bandwidth  measures  corresponding  to  a  time-limited  weighting  func¬ 
tion  are  possible  for  the  window.  Here,  we  will  consider  three  different  measures  and 
show  how  simply  they  can  be  computed  directly  from  a  set  of  discrete-time  weights 
without  having  to  compute  the  frequency  domain  window  itself. 


We  presume  a  real  and  even  set  of  weights  {wn},  such  that  =  Wn,  this 
is  different  from  the  weights  in  the  main  text,  which  were  one  sided  and  ranged 
from  iy(0)  to  w{N  —  1).  The  total  number  of  nonzero  weights  is  AT,  which  is  an  odd 
number  here.  For  time  sampling  increment  A,  the  corresponding  window  is 

W{f)  =  ^  Wn  exp  (-227r/An)  for  all  /,  (B-1) 

n 


which  has  period  1/A  in  frequency  /.  Also,  W{f)  is  even  and  real  for  all  /.  For 
future  use,  we  define  the  Nyquist  frequency  F  =  1/(2A),  which  is  halfway  to  the 
first  aliasing  lobe  at  /«  =  1/ A. 


The  first  bandwidth  measure  is  the  effective  bandwidth  fe,  defined  as  the 
single-sided  (positive  frequencies)  quantity 


(JdfW{f)j 

UfWW 


2  /  dfw(fy 


(B-2) 


However,  these  integrals  can  be  expressed  in  the  weighting  domain  very  simply  as 


I  dfW(f)  = 


-F 

F 


I  dfwu)^  = 


df  exp  {-i2'irfAn)  =  —wq  , 


exp  [— i27r/A(n  -f  m)] 


i  ^  w-n  = 


n 


(B-3) 


(B-4) 
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Substitution  in  equation  (B-2)  yields 


fe  = 


1 


n 


(B-5) 


Since  there  are  a  total  of  N  weights,  taken  at  time  increment  A,  we  can  define  the 
weighting  time  duration  as  L  =  A^A,  leading  to  the  dimensionless  product 


Lh  = 


N 

2 


n 


(B-6) 


Two  additional  “rectangular”  single-sided  bandwidth  measures  that  could 
alternatively  be  used  are  given,  by  means  of  equations  (B-3)  and  (B-4),  as 

F 

0  n 

=  (B-8) 

°  (?””j 

Observe  that  frequency  measme  /i  is  the  geometric  mean  of  the  other  two  measures; 
that  is,  fi  =  \/ /2/e)  whether  obtained  from  the  frequency  domain  or  the  time  domain 
relations. 

An  example  of  these  three  bandwidth  measures  for  the  Dolph-Chebyshev 
window,  with  total  number  of  samples  N  =  127,  is  shown  below  in  figure  B-1,  with 
the  actual  numerical  data  supplied  in  table  B-1.  The  values  of  Lfe  are  computed 
firom  equation  (B-6);  the  values  for  L/i  and  L/2  are  computed  from  the  equations 
obtained  by  substituting  L  =  NA  into  equations  (B-7)  and  (B-8),  respectively: 


r  _N  Wq 

^Jl  a  ’ 

n 

(B-9) 

N 

=  V- 

(B-10) 

(E«>n) 

B-2 


We  observe  in  figure  B-1  that  all  three  bandwidth  measures  happen  to  be  equal  when 
the  sidelobe  level  is  approximately  -23  dB.  (Note  that  all  these  ordinate  numbers  in 
figure  B-1  would  have  been  twice  as  large  if  we  had  defined  double-sided  bandwidth 
measures  in  equations  (B-2),  (B-7),  and  (B-8).) 

There  is  a  reason  for  the  upturn  in  the  L/2  curve  for  the  sidelobe  levels 
near  -20  dB.  For  deep  sidelobe  levels,  the  squared  window  W^{f)  generates  most  of 
its  contribution  to  the  integral  in  equation  (B-8)  from  the  mainlobe  region  about 
/  =  0.  However,  for  shallow  sidelobe  levels  (e.g.,  -20  dB),  where  the  mainlobe  width 
is  rather  narrow,  a  significant  contribution  to  equation  (B-8)  is  also  made  by  the 
extended  sidelobe  region.  This  causes  the  numerator  of  equation  (B-8)  to  increase, 
even  though  the  mainlobe  is  becoming  narrower.  An  alternative  bandwidth  measure, 
such  as  the  half-power  bandwidth,  would  not  have  this  upturn  for  the  shallower 
sidelobe  levels. 


Table  B-1.  The  Three  Bandwidth  Measures  Lfe,  Lfi,  and  L/2  as  a 
Function  of  the  Sidelobe  Level  for  the  127-Point  Dolph- 
Chebyshev  Window 


dB 

Lfe 

Lh 

L/2 

1^ 

warn 

Lfi 

Lf2 

-20 

0.454 

0.593 

0.774 

masm 

-22 

0.575 

-54 

lifcBia 

-24 

0.687 

0.655 

0.624 

IHI 

1.393 

1.015 

0.739 

-26 

0.784 

0.683 

0.596 

1.420 

1.033 

0.751 

-28 

0.864 

0.711 

0.584 

-60 

■WiEM 

0.931 

0.737 

0.583 

-62 

■HI 

OH 

-32 

0.988 

0.762 

0.588 

-64 

1.496 

1.086 

0.788 

mmirn 

0.786 

0.597 

-66 

Baa 

BWitel 

0.810 

0.608 

-68 

1.545 

1.119 

0.811 

-38 

1.118 

0.833 

0.620 

HI 

1.569 

1.136 

wtmm 

Bnl 

1.154 

KM 

1.592 

1.152 

-42 

1.188 

0.877 

0.647 

-74 

0.898 

0.660 

-76 

n^i 

-78 

IHI 

QH 

0.866 

-48 

1.281 

0.938 

0.687 

-80 

1.683 

1.215 

0.877 

-50 

1.310 

0.958 

0.701 

B-3 


Time-bandwidth  measure 


Peak  sidelobe  (dB) 


Figure  B-1.  Bandwidth  Measures  for  the  127-Point  Dolph-Chebyshev 
Window  as  a  Function  of  the  Sidelobe  Level  in  Decibels 
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APPENDIX  C 

KAISER-BESSEL  PROPERTIES  FOR  CONTINUOUS  WEIGHTING 


Weighting  w{t)  = 


Window  W  (/)  = 


-/o 

sin 


for  |t|  <  — ;  zero  otherwise. 


sinh  ( 


for  all  /. 


Effective  Bandwidth: 


JdfWif) 


fdfW^(f) 


[fdfW(/f 

7fdfW^(f) 


w^{0) 

2  f  dtvP{t) 


IIW 


2  J  d9  cos  (6)  Iq  {P  cos  9) 
0 


Half-Power  Bandwidth: 


W{f)  _  sinh(x)  P 
W  (0)  X  sinh  P 


,  where  x  =  \/ P^  — 


Define  fh  from  the  equation 


Wjh)  1 
W{0)  V2  ■ 


Solve 


sinh(x,,)  1  sinh(^)  rr.^  r  ^  I 

= V2~r  = V  ^ 
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First  Zero  Crossing: 


TT 


Lh 

Ml 

2 


/?2  +  TT^ 


TT^ 


V^/?2  +  TT^ 

27r 


The  three  bandwidth  measures  L/e,  L//i,  and  are  listed  as  a  function 
of  the  Kaiser-Bessel  parameter  in  table  C-1.  These  data  are  also  plotted  in  figure 
C-1.  Table  C-1  and  figure  C-1  convey  the  same  information  for  the  Kaiser-Bessel 
weighting  as  figure  B-1  and  table  B-1  do  for  the  Dolph-Chebyshev  weighting. 


Table  C-1.  The  Three  Bandwidth  Measures  Lfe,  Lfh,  and  \Lfz  as  a 
Function  of  the  Window  Parameter  j3  for  the  Kaiser-Bessel 
Window 


p 

Lfe 

Lfh 

\Lf, 

2 

0.7544 

0.4966 

0.5927 

3 

0.9398 

0.5467 

0.6914 

4 

1.1006 

0.5997 

0.8095 

5 

1.2396 

0.6516 

0.9398 

6 

1.3634 

0.7011 

1.0779 

7 

1.4763 

0.7481 

1.2211 

8 

1.5809 

0.7928 

1.3679 

9 

1.6789 

0.8352 

1.5172 

10 

1.7713 

0.8758 

1.6682 
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APPENDIX  D 

MATLAB  CODE  TO  COMPUTE  EDF 


This  appendix  lists  the  MATLAB  code  that  can  be  used  to  compute  the 
EDF  for  a  given  windowing  function.  This  code  is  available  from  the  authors  upon 
request. 

function  [Edf,  Fo]  =  edf(W,  N,  T,  Pm) 


y,  edf.m  -  evaluates  the  expected  degrees  of  freedom  for  a  given 
y.  windowing  function  for  various  values  of  the  percentage 

y.  overlap.  The  formeil  parameters  are  defined  as  follows. 

y. 


y. 

y. 

y. 

y. 

y. 


W  -  temporal  weights  (N  x  1) 

N  -  size  of  the  FFT 

T  -  total  number  of  data  points 

Pm  -  maximum  number  of  pieces  to  be  used 


%  Edf,  a  Pm  X  1  vector,  is  returned  to  the  calling  context, 

y«  as  well  as  the  corresponding  values  of  the  fractional  overlap 

y,  in  the  Px  x  1  vector  Fo. 

y. 

%  Edf  -  expected  degrees  of  freedom  (Pm  x  1) 

%  Fo  -  fractional  overlap  (Pm  x  1) 

y.  $Id:  compute.edf  .tex,v  1.2  1996/09/23  07:23:48  hall  Exp  hall  $ 


Edf  =  zeros (Pm,  1); 

Fo  =  zeros (Pm,  1); 

S  =  zeros (Pm,  1); 

Phi2  =  zeros (N,  1); 

N1  =  N  -  1; 

PhiO  =  W’  *  W;  %  dot  product  of  the  window  W 

P2  =  PhiO  *  PhiO; 

Phi2(l)  =  1.0; 

for  Ks  =  1:N1 
A  =  0.0; 
for  Ns  =  Ks:Nl 
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A  =  A  +  W(Ns  +  1)  *  W(Ns  -  Ks  +  1); 

end 

Phi2(Ks  +  1)  =  A  *  A  /  P2; 

end 

Edf(l)  =  2.0; 

y,  P  is  the  number  of  pieces  of  data  being  used.  All  of  the  values 

y.  of  EDF  for  P  less  than  (T/N)  are  the  same  because  there  is  no  overlap. 

for  P  =  (T/N): Pm 

Sc  =  (T  -  N)  /  (P  -  1) ; 
for  Ps  =  1:P 

S(Ps)  =  round((Ps  -  1)  *  Sc); 

end 

A  =  0.0; 

for  Ps  =  1 : (P-1) 

Sp  =  S(Ps); 
for  Qs=(Ps+l):P 

Ns  =  S(Qs)  -  Sp; 
if  (Ns  <  Nl) 

A  =  A  +  Phi2(Ns  +1); 

end 

end 

end 

Edf(P)  =  2  *  P  *  P  /  (P  +  2  *  A);  */.  EDF  from  equation  (13) 

Fo(P)  =  max(0,  1  -  Sc  /  N) ;  */.  FO  from  equation  (12) 

end 

y,  end  of  edf.m 
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